Dynamic instability of speckle patterns 
in nonlinear random media 
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Linear stability analysis of speckle pattern resulting from multiple, diffuse scattering of coherent 
light waves in random media with intensity-dependent refractive index (noninstantaneous Kerr non- 
linearity) is performed. The speckle pattern is shown to become unstable with respect to dynamic 
perturbations within a certain frequency band, provided that nonlinearity exceeds some frequency- 
dependent threshold. Although the absolute instability threshold is independent of the response 
time of nonlinearity, the latter significantly affects speckle dynamics (in particular, its spectral con- 
tent) beyond the threshold. Our results suggest that speckle dynamics becomes chaotic immediately 
beyond the threshold. 
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I. INTRODUCTION 

Recent years have been marked by a considerable 
progress in understanding of multiple scattering of light 
in random media in the so-called mesoscopic regime, 
when the coherence length Lcoh of the incident light ex- 
ceeds both the size of the disordered sample L and the 
absorption length and Lcoh, La ^ I, where I 
is the mean free patbi*2i^. In the mesoscopic regime, 
interference of multiply scattered waves gives rise to a 
complicated and seemingly random spatial distribution of 
light intensity called 'speckle pattern' or just 'speckle' for 
brevity. Speckles exhibit a number of interesting prop- 
erties, such as short- ^ and long-ranga^iSiLSiSiifliiLiSii^ 
correlations (and related memory effecfei^ and 'univer- 
sal conductance fluctuations'—), and enhanced sensi- 
tivity to small displacements of scatterersiSiiLi^. As 
a practical matter, multiply scattered light finds ap- 
plications in optical characterization of turbid mate- 
rials (e.g., diffusing- wave spectroscopjiiSiii and optical 
microrheology^-), medical diagnostics^^, design of ran- 
dom laser sourcesSiiS^, and cryptographji^^. 

Most of the available work on multiple light scatter- 
ing and speckles has concentrated on linear random me- 
dia, where the polarization of the medium depends lin- 
early on the electric field of the light wave. This regime 
is most often realized in practice since optical nonlin- 
earities are known to be (generally) weak24S^. At the 
same time, nonlinear optical systems are known to ex- 
hibit much more complex and reach behavior than their 
linear counterpartsS^i2SiSS42L2&, and it would be there- 
fore important for both fundamental research and ap- 
plications to advance our understanding of optical phe- 
nomena in nonlinear random media. Up to now, the 
interplay of randomness and nonlinearity has been stud- 
ied mainly in one-dimensional systems, where numerical 
methods proved to be rather efficient (see, e.g., Refs. 
andls^ for representative bibliography). Because waves 
are known to behave quite differently in systems of dif- 
ferent dimensionalities, one-dimensional results are not 



readily generalizable to the three-dimensional case. In 
three-dimensional nonlinear random media, certain at- 
tention has been paid to the behavior of average trans- 
mission characteristics in the contexts of the so-called 
self-transparency effecl— and optical power limitation—, 
whereas available studies oi fluctuations of multiply scat- 
tered waves are rather limited^'^i^'^i'^^i^^i^'''i^^ . 

Recently, a new phenomenon — dynamic instabil- 
ity of speckle patterns in random media with intensity- 
dependent refractive index (Kerr nonlinearity) — has 
been theoretically predicted^S*^. Theoretical study of 
speckle dynamics beyond the instability threshold has 
been initiated in Ref. using a newly developed dy- 
namic Langevin approach. It has been found that speckle 
pattern is likely to exhibit a transition from a station- 
ary state to the chaotic dynamics upon increasing the 
strength of nonlinearity. The analysis of Ref. suffers, 
however, from numerous limitations, making it hardly 
applicable to real experimental situations. Namely, the 
nonlinearity is assumed to be instantaneous (nonlinear- 
ity response time tnl = 0), only the limit of very large 
sample size {L/l 3> where k = 27r/A and A is the 
wavelength of light) is treated, and the typical time scale 
T of speckle dynamics beyond the instability threshold 
is assumed to be quite large: r ^ {£/c){k£)'^, where 
c is the speed of light in the average medium. Mean- 
while, noninstantaneous nature of nonlinearity is known 
to play an important role in nonlinear dynamics of op- 
tical systema^^i-^3,44^45,46,47^ ^ijg condition L/£ > fc£ may 

be difficult to fulfill in practice, and imposing severe con- 
ditions on T in the beginning of analysis is undesirable 
since r is one of the quantities to be found as a result. 

In the present paper, we apply the dynamic Langevin 
approach introduced in Ref. 41 to perform a detailed 
study of the onset of speckle pattern instability and 
speckle dynamics beyond the instability threshold, re- 
laxing the above limitations''®. We show that the insta- 
bility threshold is actually lower than found in Ref. 
coinciding with the result obtained in Refs. [s^ and l40l 
in a completely different way. Dynamics of speckle pat- 
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tern beyond the instability threshold appears to be very 
sensitive to relations between k£ and L/£, and t^l and 
Td, where Td = L^/D is the time required for multiply 
scattered light to propagate through a random sample 
of size L, and D is the photon diffusion constant. Al- 
though we expect the speckle dynamics to be chaotic be- 
yond the threshold independent of particular values of 
above parameters, the spectral content of this dynamics 
and, in particular, the maximal excited frequency and 
Lyapunov exponents, are found to be sensitive to the ra- 
tios k£/(L/£) and tnl/^d- Finally, we comment on the 
accuracy of our theory by identifying physical processes 
neglected in the calculation and discuss the conditions of 
possible experimental observation of the instability phe- 
nomenon for light in nonlinear random media. 

II. THEORETICAL MODEL 

We consider multiple scattering of monochromatic (fre- 
quency Lu) light in a non-absorbing sample of random 
medium (typical size L) with weak disorder (k£ ^ 1) 
and intensity-dependent dielectric constant. For simplic- 
ity, and taking into account that multiple scattering de- 
polarizes the light, we apply scalar approximation from 
the beginning. The problem then reduces to the following 
scalar nonlinear wave equation: 

- ^£ [1 + '^^W + ^^NL(r, t)]| Eir, t) 

= J{r,t), (1) 

where E{r, t) is the wave amplitude, J(r, t) — 
Jo(r) exp(— icji) is a monochromatic source term, Se{r) 
is the fluctuating part of the linear dielectric constant 
normalized by its average value (e), and the nonlinear 
part of the dielectric constant AeNL(r,^) (normalized by 
(e) as well) depends on the intensity I{r,t) = \E{r,t)\^ 
of the wave. For de{r) we adopt the Gaussian white- 
noise model: (&(r)fc(ri)) = in/ {k'^£)5{r — ri). It will 
be clear from the following that the specific model used 
for Se{r) should not affect the main results of the paper 
significantly because they originate from the long-range 
diffusion of light in the random sample and not from the 
short-range properties of dielectric constant fluctuations. 

If the random medium were linear [AeNL(r,t) = 0], 
the time-dependence of E in Eq. would be har- 
monic [inasmuch as the source term J(r, t) is harmonic] : 
E{r,t) = Eo{r) exp{~iujt), and hence / = |£'o(r)| would 
be independent of t. Due to nonlinearity, however, the 
optical system under study can exhibit instabilities re- 
sulting in spontaneous (periodic or chaotic) dynamics of 
/, despite the fact that the incident radiation [Jo(r)], 
nonlinear sample [(5£(r)], and conditions of light prop- 
agation remain unchanged (i.e., are time-independent). 
Instability phenomena are widely encountered in non- 
linear optic£^^^^i^^i4'^i^^^-^i-^5,46,47_ Allowing for sponta- 
neous dynamics in our random optical system (but in no 



way imposing its existence a priori), we keep the time de- 
pendence of / and Aenl in what follows, assuming, how- 
ever, that this dependence is relatively slow, i.e. that 
neither / nor Aenl do not change significantly on the 
time scale of the mean free time i/c. It will be seen from 
the following that spontaneous wave dynamics in non- 
linear random media is indeed possible for a sufficiently 
strong nonlinearity and that a regime exists when the 
characteristic time scale of this dynamics is much larger 
than £/c. 

Due to the randomness of 6e{r), E{r,t) and I{r,t) 
are random functions of r as well (speckle pattern). We 
will be therefore interested in their statistical moments 
rather than in their specific r- and i-dependences for 
some particular realization of Ss(r), an approach com- 
mon in statistical 0pticcii'^iVi^'^i^i«i"ii4ii5ii6,i7,i8,i9^ 

particular, we will analyze the correlation function of in- 
tensity fluctuations {6I{r,t)SI{ri,ti)), where 61{r,t) — 
/(r, i) — (/(r)), angular brackets denote averaging over 
realizations of disorder fe(r), and (/(r)) is the average 
intensity that we assume time- independent. In linear me- 
dia, spatial variations of are known to be smooth; 
(/(r)) changes significantly on the scale of the sample 
size L but not on the scale of the mean free path £. In 
the following we assume that this smoothness of (/(r)) 
is preserved in a nonlinear medium. To further sim- 
plify the analysis, we assume that the nonlinearity is 
weak enough to cause no significant modification of the 
mean free path £. It can be showr.'^^ that this requires 
AeNL(r, i)|^ fc£ <JC 1. On the one hand, if this condition 
is violated, our analysis still remains at least qualitatively 
valid, provided that the mean free path £ is recalculated 
with account for nonlinearity of the medium. On the 
other hand, the condition \Ae-M-L{j^,t)f k£ ^ 1 does not 
really limit the validity of our main results because, as 
will be seen from the following, instability of speckle pat- 
tern is expected to develop for much weaker nonlinearities 
than those required to violate the above condition. 

It is known that in a linear medium, correlation func- 
tion of intensity fluctuations {SI{r, t)5I{ri,ti)) comprises 
short- and long-range contributions. The short-range (in 
space) contribution (often denoted by Ci), 

{SI{r,t)SI{r,,h)),^,,, = \{E{r,t)E{r,h))\^ 

X /(r-rl)^ (2) 

where /(r) = [sin(fcr)/(A;r)] exp[— r/(2^)], is large for 
r — ri| < £ but decreases to zero rapidly at larger sep- 
arations between r and ri— . The-long range contribu- 
tion is much weaker [oc {k£)~^] but persists through the 
whole random sample (i.e. even for jr — ri| ~ L )^i^/ . 
It can be theoretically calculated by using the Langevin 
approach^'^ that amounts to writing the Langevin equa- 
tion for SI{r, t), 

^5I{v,t)-DV^5I{v,t) - - V •jext(r,i), (3) 
at 

where jcxt(r,t) are random Langevin currents 
that have zero mean and correlation functions 
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Here i,j = x,y,z and j^^^ denotes the corresponding 
Cartesian component of joxt- Short-range correla- 
tion function of intensity fluctuations and correlation 
function of Langevin currents can be schematically 
represented by the diagram (a) of Fig. ^ Exact func- 
tional dependence of the long-range correlation function 
((5/(r, i)^/(ri, on the separation Ar = r — ri 

between the points r and ri, as well as on the time 
variables, depends on the geometry of disordered sample 
and source of waves. In the infinite random medium 
with (/(r)) « /q and for t = ti one finds^ 

(W(r,t)^/(ri,t)),„„^^^. (4) 

The spatial long-range correlation function of intensity 
fluctuations (often designated as C2) has been mea- 
sured experimentally for microwaves in a quasi-one- 
dimensional disordered waveguideSiiS. In that case, the 
magnitude of correlation is of the order of 1/g (where 
g 3> 1 is the dimensionless conductance) and not of the 
order of {M)^'^ [as in Eq. Q], but it exhibits a similar de- 
pendence on Ar as in three dimensions: C2 oc 1/Ar. Spa- 
tial correlation function of intensity fluctuations is known 
to have infinite-range components often designated as C3 
(see, e. g., R ef.Eland references therein) and Co (see Refs. 
Ha and Hal . We do not take these latter terms into ac- 
count in what follows because C3 oc {kt}~^ -C C2 and Co 
takes considerable values only for the case of spatially- 
localized ('point-like') sourcei^, whereas we assume an 
extended source of light. 

In linear media, Langevin currents joxt(r,t) are very 
sensitive to changes of disorder Seir)^. Generaliz- 
ing this property, one concludes that in a nonlinear 
medium joxt(r,i) should be sensitive to changes of 
fe(r) -I- A£NL(r,t)^, and, more precisely, to changes 
of AeNL(i", since 8e{v) is fixed (and hence does not 
change) in our case. Modification Ajoxt(r,i) of Langevin 
currents due to an infinitesimal change A[AeNL(r, i)] of 
the nonlinear part of dielectric constant is^ 

Ajcxt(r,t) - j d?v' [ dt'q(r,r',t-t') 

Jv J -00 

X A[AeNL(r',t')], (5) 

where the spatial integral is over the volume V of the 
sample, we neglect the terms of the second and higher or- 
ders in AeNL(r, t), and q(r, r',i— i') are random, sample- 
specific response functions describing the sensitivity of 
Langevin currents jcxt at r at time t to small changes 
of dielectric constant e = 1 + 6e + Aenl at r' at time 
t'. (In general case, A[A£NL(r', i')] Eq. (|3J) should 
be replaced by A£(r',t'), but we assume that the lin- 
ear dielectric constant 1 -I- Se{r) is fixed.) The physical 
meaning of q can be understood by considering a brief 
'fiash' of e, Ae(r,t) oc 5{r — Tq) S{t — to), that occurs at 
time to and is localized at Tq. According to Eq. |SJ|, for 




FIG. 1: (a) Diagrammatic representation of the correlation 
function (^Jcxt(r, i)jcil(ri, ii)^ of Langevin currents jcxt (r, t). 
(b)-(d) Diagrams contributing to the correlation function 
g*-'-* (r, r ', At)?'-*' (ri, rj", Afi)^ of random response functions 

q(r,r',At). (e) and (f) are examples of diagrams neglected 
in our analysis. Solid lines with right (left) arrows denote re- 
tarded (advanced) Green's functions of the linear disordered 
wave equation. Dashed lines denote scattering of the two con- 
nected wave fields on the same heterogeneity. Wavy lines are 
fcg vertices. The diagrams (b)-(f) are obtained by functional 
differentiation of the diagram (a) with respect to the dielectric 
constant of random medium. 

t > to such a fiash will modify the Langevin current at 
some distant point r by Ajoxt(r, t) oc q(r, ro,t — to). To 
account for changes of e throughout the whole disordered 
sample and at all times before the time t, one has to sum 
their associated Ajcxt, which results in integration over 
space and time and leads to Eq. (Q. 

Obviously, the mean value of q is zero, while the cor- 
relation function of its Cartesian components is given by 
a sum of diagrams (b)-(d) of Fig. Q 

(g«(r,r',At)gW(ri,r{,Ati)) 
^3T:D^{c^/£)S,jd{r-ri) 

X {(/(r')) G(r', r[;At - Aii)G(r{, r; Aii) (/(r)) 

+ (/(rO) G(r{, r'; At, - At)G(r', r; At) (/(r)) 

- (/(r'))G(r',r;At) {I{r[)) G{r[,r; At,) 

+ ei/{inD)5{At - A<i)/2(r' - y[) (/(r)) 

X [(/(r'))G(r',r;At) + (/(rO)G(r{,r;At)]}. (6) 

Each term in Eq. ^ and a corresponding diagram of Fig. 
n can be given a simple physical interpretation. Con- 
sider, for example, the first term, (/(r ')) G(r ', r{; At — 
Ati)G(r{,r; A^i) (/(r)), given by the diagram (b) of Fig. 
n This term contains two wave fields E and two complex 
conjugate fields E* that enter the disordered medium 
from the left. The first pair of E and E* , shown by two 
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lower solid lines in Fig. QJb), propagate together to the 
vicinity of points r and ri, experiencing scatterings at 
exactly the same scatterers (this is denoted by dashed 
horizontal lines). After the last scattering event, these 
two fields do not follow the same path anymore: E goes 
to ri and E* goes to r. The second pair of E and E*, 
depicted by two upper solid lines in Fig. , propagate 
together to some intermediate point r' where E experi- 
ences scattering on a fluctuation A[AeNL(r', ^ ')] of the 
nonlinear part of dielectric constant. The two wave fields 
then continue together to some other point r{ , where E* 
is scattered by A[AeNL(i"i, ^i)]- Finally, the two wave 
fields make the remaining way from rl to the vicinity of 
points r and ri together and then separate for E to go 
to r and E* — to ri. Similar interpretation can also 
be given to other terms of Eq. © and to the diagrams 
(c)-(f) of Fig. n We note that the last term in Eq. © 
[the term containing the factor /^(r' — r[)], given by 
the diagram (d) of Fig. has been neglected in Ref. l4ll 
This is only justified in the limit of extremely large sam- 
ple size and for very slow speckle dynamics as will be 
shown below. Also, it should be noted that there exist 
other diagrams that contribute to the correlation func- 
tion (gW(r,r', At)g(j)(ri,r{, Ail)) but are neglected in 
our analysis. We show examples of such diagrams in pan- 
els (e) and (f) of Fig. ^ These two particular diagrams 
are neglected in Eq. © because their contribution to our 
principal result, Eq. (|12() . is a factor {£/L)'^ <C 1 smaller 
than the contribution of diagrams (b)-(d). In addition to 
neglecting non-leading diagrams in Eq. © , we note that 
the diagrams (b) and (c) will enter further analysis with 
the assumption |r' — r{| > £. Taking into account the 
opposite possibility (|r' — r{| < £) is possible but would 
modify our final Eq. p2(l only by additional terms that 
are a factor £/L ^ 1 or a factor {k£)~^ ^ 1 smaller than 
the terms due to |r' — r{| > £. 

A dynamic equation for joxt(r, t) can be obtained from 
Eq. lO by assuming that the changes of Langevin cur- 
rents and nonlinear dielectric constant entering into that 
equation take place during some time interval St. We 
then divide both sides of the equation by St and take the 
limit oi St ^ 0. This yields 

d [ r°° 

— jext(r,t) = y^^'^y^ dAtq(r,r',Ai) 

X ^AeNL(r',t- Ai). (7) 
at 

To describe the temporal evolution of AeNL(r, t), we em- 
ploy a Debye-type relaxation equation, 

d 

tnltt: AeNL(r, i) = "AeNL(r, t) + 2n2«(r, t), (8) 
at 

which is commonly used in nonlinear optics to model 
the noninstantaneous medium response^^'^S'^Sii^iii. Here 
the constant n2 measures the strength of nonlinearity 
and the relaxation time tnl characterizes the speed of 
nonlinear response of the medium. Note that we have 



subtracted from AeNL(r,i) in Eq. ^ its average value 
2n2 (/(r)) ^ 1 that only leads to a small modification 
of the average refractive index of random medium and 
is of no significance in the context of the present study. 
Equations lO and (jS]) form a self-consistent system. 
Equation Q governs the spatio-temporal evolution of 
the intensity fluctuations 5/(r, t) due to the Langevin 
currents jcxt(r, i), Eq. {Tj) describes the distributed feed- 
back mechanism, leading to variations of jext(r, i) de- 
pending on the changes of AeNL(r,t), and Eq. © re- 
lates A£NL(r,*) to SI{T,t) = /(r,t) - (/(r)). Note that 
Eq. |(7J) is a linearized equation: only the terms linear 
in AeNL(r,i) 1 are kept. As we show below, the lin- 
earized nature of Eq. Q may result in the exponential 
growth of its solution with time, and in this sense Eqs. 
(0) , O , and © are analogous to the equations of linear 
stability analysis commonly used to study the stability 
of nonlinear systems^^-^^. Hence, although Eqs. (O, iQ, 
and ||HJ) will allow us to study the stability of the speckle 
pattern and the typical time scales of spontaneous inten- 
sity fluctuations beyond the instability threshold, they 
cannot be used to determine the amplitude of these fluc- 
tuations. 



III. LINEAR STABILITY ANALYSIS 

Consider now an infinitesimal periodic excitation of the 
stationary speckle pattern: (5/(r, i) = SI{y,v) exp(wi), 
where v = ^1 — lA. Such an excitation can be either 
damped or amplified, depending on the sign of the Lya- 
punov exponent A. The value of A is determined by 
two competing processes: on the one hand, diffusion 
tends to smear the excitation out, while on the other 
hand, the distributed feedback sustains its existence. The 
mathematical description of this competition is provided 
by Eqs. (0), and ijHl that after the substitution of 
SI{r,t) = SI{r,v) explii/t) [and similarly for joxt(r, i) and 
A£NL(r,i)] become 

tvSIir, v) - DV^SI{v, v) = -V- jext(r, v), (9) 

ivicAr,iy) ^iv d^r' / dAi q(r, r ', At) 
Jv Jo 
xAeNL(r',t^)exp(-wAt), (10) 

ij/rNLA£NL(i", 1^) = -A£nl(i", i^) + 2n2(5/(r, v). (11) 

If = 0, Eq. Q yields the stationary, time- independent 
component of the speckle pattern, Eq. I|1U|I is trivial, and 
Eq. (inj deter mines the time-independent part of the 
nonlinear dielectric constant (more precisely, its devia- 
tion from the average). This case is not of interest for us. 
From here on, we will consider the case of ^ 0, because 
we are interested in dynamics of speckle pattern. We pro- 
ceed by expressing Ae^i^{r,u) through SI{r,v) by using 
Eq. Hll|l and the correlation function {SI{r, i')SI*{ri, v)) 

through ^j^^^t(r,j/)ji^i*(ri,;/)^ by using Eq. ©. The 

latter correlation function can, in its turn, be expressed 
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through {SI{r,v)SI*{ri,v)) by using Eq. (fTUIl . The two 

resulting equations relating ^j^^j(r, j/)j^:^|.*(ri, j/)^ and 

{6I{r,i/)5I*{ri,i')) have to be consistent, which yields 
an equation for allowed values of H and A. The calcula- 
tion can be performed along the same lines as in Ref. ^3 
and yields 

l + k£/(L/£) 



h{nTD,ATD) + kl/(L/l)g{KTD) 
X i/(r!rNL,ArNL), (12) 



where An = 7i2(/(r)) is a typical value of the nonlin- 
ear modification of refractive index (with bar denoting 
spatial average) and the dimcnsionless function h is de- 
fined in Ref. liJ [Eq. (AlO) of Ref. ^ with R = 0]. In- 
stead of reproducing here a (rather lengthy) expression 
for /i(x, y), we only note that /i(0, 0) = 1, 1 — h{x^ 0) cx 
for a; ^ 1 and h{x,Q) oc a;~^/^ for a; ^ 1. Functions H 
and g in Eq. (|12|l are defined as -ff (x, y) = x"^ + {\ + y)^ 
and 



giy) 



1 



(2-3e-i/2)2y+l 

exp(^-^v/2yTT^ 



1 



(13) 



Finally, in Eq. ^ we have introduced a bifurcation pa- 
rameter 



L 



(14) 



As will be seen from the following, this parameter con- 
trols the stability of speckle pattern with respect to weak 
perturbations. 

Different terms in Eq. H12I) have precise physical mean- 
ing. Functions h and g describe the distributed feedback 
mechanisms due to long-range [h) and short-range [g) 
correlations of intensity fluctuations in random medium. 
As follows from Eq. H12|) . the relative role of these 
two feedback mechanisms is controlled by the parame- 
ter kij {L/t). The fact that the nonlinear response of the 
medium is noninstantaneous is described by the function 
H . H = 1 for tnl = and H increases monotonically 
with both Vl and A for tnl > 0. 



IV. INSTANTANEOUS NONLINEAR 
RESPONSE 

We first illustrate the main result of the linear stabil- 
ity analysis performed in the previous section — Eq. \12\ 
— on the most simple case of instantaneous nonlinear- 
ity, assuming tnl = 0. In this case, H (firNL, Atnl) = 1 
in Eq. (|12|l . For every value of p and a given Q. we can 
then solve Eq. \i2\ for A. If A < 0, a weak harmonic 
(in time) excitation of speckle pattern at the chosen fre- 
quency Vl is exponentially suppressed and hence disap- 
pears after a time ^ 1/ |A|. If, in contrast, A > our 




FIG. 2: Stability diagram of speckle pattern in a nonlinear 
disordered medium. Physically realizable combinations of ex- 
citation frequency n, Lyapunov exponent A and bifurcation 
parameter p belong to the surface shown in the figure. Thick 
line corresponds to A = and splits the surface into two parts. 
Instability (A > 0) can only develop in the upper part of the 
surface and hence requires p > 1. In the lower part of the 
surface, p < 1 and A < for all fi. 



analysis predicts exponential amplification of the excita- 
tion and hence instability of speckle pattern with respect 
to weak perturbations at frequency 17. It follows from 
Eq. that A < for all frequencies fi, if p < 1. We 
illustrate this in Fig. |2 where we show a surface plot of 
p versus VI and A for a particular case of k£/{L/£) = 1. 
A thick line in Fig. |21 shows a crossection of the surface 
by a plane A = and splits the surface into two parts. 
All points with A > lie in the upper part of the sur- 
face and correspond to p > 1. At higher f2, larger p is 
required to get into the upper part of the surface. All 
points corresponding to p < 1 belong to the lower part 
of the surface and correspond to A < 0. When p exceeds 
1, A becomes positive in a band of frequencies (0, fimax) 
and speckle pattern is hence expected to exhibit insta- 
bility with respect to low-frequency excitations. Thus, 
p = 1 is the absolute instability threshold. This conclu- 
sion agrees with the result of Refs. !33 and '2(3 (where the 
same expression for the absolute instability threshold has 
been found by using a different approach) and coincides 
with the result of Ref. 41 in the limit oi L / £ ^ k£ con- 
sidered in that paper. Equation (|12l) allows us to study 
dynamics of speckle pattern at p > 1 and arbitrary L/£, 
H> 1. 

Frequency-dependent instability threshold (i.e. the 
value of p at which A = for a given H.) is shown in Fig. 
13 for several values of k£/{L/£). If k£/{L/£) < 1 and 
n fli = T^^[{L/£)/Jk£)f, which corresponds to the 
limit considered in Ref. 41;, the denominator of Eq. (|12|) is 
dominated by the first term h{QT]:), ATo) and hence the 
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FIG. 3: Frequency dependence of instability threshold in a 
nonlinear disordered medium with instantaneous Kerr non- 
linearity (rNL/ro ~ 0) for k(./{L/£) — > (upper solid curve) 
and k£/{L/£) — 0.1, 1, and 10 (solid curves from top to bot- 
tom). Inclined dashed lines show scaling laws D,^ and fi^''^. 
Horizontal dashed lines show asymptotic values of p for large 

n. 



FIG. 4: Lyapunov exponent A versus excitation frequency 
fl (both in units of Tj,"^ = D'^/L) in a nonlinear disordered 
medium with instantaneous Kerr nonlinearity (tnl/Td ~ 0) 
slightly above the absolute instability threshold (p — 1.1) for 
M/{L/£) (lower solid curve) and M/{L/e) = 0.1, 1, and 
10 (solid curves from bottom to top). Horizontal dashed line 
shows A = 0. 



threshold value of p obeys scaling laws p ~ 1 cx {QTy^)^ 
for n < To ^ and p oc (f^TD)^/^ for TjJ^ < 1] < Vti [see 
the curves corresponding to kt/{L/t) — 0, 0.1 and in- 
clined dashed hues in Fig. |3] . This regime is governed by 
the presence of long-range correlations of intensity fluc- 
tuations [see Eq. and is mathematically described by 
the diagrams (b) and (c) of Fig. ^ The upper hmiting fre- 
quency of the band inside which the speckle pattern is un- 
stable (instabiUty band) scales as fimax oc T^^{p — 1)^/^ 
for < p — 1 -C 1 and as ilmax cx; T^^p^ for p — 1 > 1. 
Note that Q.i = i[(L/^)/(H)]2 - {c/£)/{Mf exceeds 
^ for Lli> M but is still a factor (fc£)'^ < 1 smaller 
than the inverse mean free time c/ £. Therefore, our anal- 
ysis applies for > ili as well. In this latter case, 
the second term, k£/{L/£)g{ATu), in the denominator 
of Eq. H12|) becomes larger than the first one and we find 
p ~ 1 -I- {L/£)/{k£) (see dashed horizontal lines in Fig. 

. As one can see from Fig. the curves corresponding 
to k£/{L/£) = 1 and 10 indeed tend to this asymptotic 
value for > Sli. If p exceeds 1 + {L/£)/{k£), A > for 
all ft and speckle pattern becomes unstable with respect 
to periodic excitations at all frequencies (in other words, 
f^max — *■ oo). It is worthwhile to note that this wide-band 
instability is due to strong short-range correlations be- 
tween intensity fiuctuations in random media [see Eq. |(5J] 
and mathematically results from the diagram (d) of Fig. 
n The short-range correlations appear to dominate the 
development of speckle pattern instability in relatively 
small random samples {L/£ < k£, but still L/£ 3> 1) 
and for high frequencies {Q > Qi), while in large sam- 
ples {L/£ ^ k£) and for low frequencies (fi < fii) the 
principal role is played by the long-range correlations. 



Lyapunov exponent A as a function of frequency f2 
can also be obtained from Eq. H12|) . For a system slightly 
above the absolute instability threshold {p — 1.1), we 
plot corresponding results in Fig. ^ for the same values 
of k£/{L/£) as in Fig.O In all cases, the maximal value 
of Lyapunov exponent Amax is reached for ^ and A 
decreases monotonically with f2. 



V. NONINSTANTANEOUS NONLINEAR 
RESPONSE 

The case of instantaneous nonlinear response, tnl = 0, 
considered in the previous section, is, obviously, not very 
realistic. In this section we extend the analysis of the 
previous section to arbitrary tnl- For concreteness, we 
fix the ratio k£/ {L/£) to a realistic value of 1 and plot the 
frequency-dependent instability threshold resulting from 
Eq. H12|) in Fig. |S1 It is readily seen from the figure that 
the noninstantaneous nature of nonlinearity does not af- 
fect the absolute instability threshold p — I but modifies 
significantly its frequency dependence. We have now sev- 
eral typical times (frequencies) that enter into play: the 
time Td needed for a wave to diffuse through the random 
sample, the time tnl of nonlinear response, and the time 
n^'^, the inverse of the frequency at which the instability 
threshold saturates at a constant level in the case of in- 
stantaneous nonlinearity (corresponding typical frequen- 
cies are TjJ^, and fli, respectively). In the particular 
case of k£/{L/£) = 1 that we consider, fii = TjJ^. 

Let us first assume that tnl ^ • For Q <^ t^^^I , the 
behavior of instability threshold is then similar to that in 
the case of instantaneous nonlinearity (see Fig. O and the 
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FIG. 5: Frequency dependence of instability threshold in 
a nonlinear disordered medium with noninstantaneous Kerr 
nonlinearity (tnl /Td = 1, 10"^ and 10"^ solid curves from 
top to bottom) compared with the case of instantaneous non- 
linearity (tnl/Td = 0, lower solid curve). kl/{L/t) — 1 for 
all curves. Inclined dashed line shows scaling law Q,^ , horizon- 
tal dashed line shows asymptotic value of p for large Q and 
tnl/Td = 0. 
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FIG. 6: Lyapunov exponent A versus excitation frequency 
Q, (both in units of TjJ^ = D^/L) in a nonlinear disordered 
medium with noninstantaneous Kerr nonlinearity (tnl/Td = 
1, 10"^ and 10~^ solid curves from bottom to top), slightly 
above the absolute instability threshold (p = 1.1), compared 
to the case of instantaneous nonlinearity (tnl/Td ~ 0, upper 
solid curve). k£/{L/£) = 1 for all curves. Horizontal dashed 
line shows A = 0. 



lower solid curve of Fig. [Sj. First, the threshold value of 
p—1 scales as {V,Ty))^ for Q < ^, then the curve starts 
to saturate at a constant level for > i7i [if L/£ > kl 
and T]^^ < there is a region T^^ < 17 < fii, where 
p oc {^TdY/'^]. However, in contrast to the case of instan- 
taneous nonlinearity, this saturation is not complete and 
is followed by a parabolic growth of p with for > Tj^^ . 
Physically, this is due to the fact that the nonlinear part 
of dielectric constant Aenl cannot follow such fast dy- 
namics and the amplitude of its oscillations is strongly 
suppressed for > t^I^. Mathematically, the growth 
of the threshold value of p with for > t^^^ origi- 
nates from the term (firNL, Atnl) in Eq. H12|) . This 
behavior is illustrated by the solid curves corresponding 
to tnl/Td = 10-3 and 10"^ in Fig.EJ 

If the response time of nonlinearity tnl is of the or- 
der or shorter than Td (see the upper solid curve in Fig. 
ISJ , the frequency dependence of instability threshold be- 
comes dominated by the noninstantaneous nature of non- 
linearity at all frequencies. We find p — 1 (x (JItnl)^ as 
illustrated in Fig. |Slfor tnl/Td = 1. 

The frequency dependence of Lyapunov exponent A 
found by solving Eq. (|12|) slightly above the absolute in- 
stability threshold (p = 1.1) and for kl/{L/t) = 1 is 
shown in Fig. for several values of tnl/Td. Qualita- 
tively, behavior of A is similar to that found in the case 
of instantaneous nonlinearity (see Fig. ^ and the upper 
solid curve of Fig.lSJ, but A decreases with Q. faster than 
for tnl/2d = 0. Together with a faster growth of in- 
stability threshold (Fig.[SJ), this narrows the band of fre- 
quencies (0, r^max) where excitations are subject to insta- 



bility. Instead of becoming infinite for large p as in the 
case of instantaneous nonlinearity (tnl = 0), rJmax scales 
as t'Ip^/'^ for tnl/Td > 1. 

VI. PATH-INTEGRAL PICTURE OF SPECKLE 
PATTERN INSTABILITY 

It is worth noting that most of the peculiarities of the 
above analysis can be qualitatively understood in the 
path-integral picture of wave propagation, especially in 
the case of instantaneous nonlinearity that we now con- 
sider in more detail. 

Consider a partial wave traveling along a typical dif- 
fuse path of length s ^ L"^ / 1 va. & nonlinear disordered 
medium. Assume that the path length s can vary slowly 
with time due to, e.g., an infinitely slow Brownian mo- 
tion of scatterers. A diffuse wave following the same 
path a time At later will then acquire a different phase. 
The phase difference between the two waves, A0, can be 
splited in a sum of two contributions: a 'linear' contri- 
bution A0L, resulting directly from the change of path 
length s and independent of the nonlinear properties of 
the medium, and a 'nonlinear' contribution A(/)nl due to 
the nonlinearity of the medium: 

A0NL = kn2 f dsiAI{r), (15) 
Jo 

where A/(r) denotes the change of intensity during the 
time At, we take into account the fact that the nonlinear 
part of the refractive index is simply n2/(r) in the consid- 
ered case of instantaneous nonlinearity, and the integral 
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is along the wave path. In the Umit of weak nonUnearity, 
it can be show n"^"-"^^. that the mean values of both A0l 
and A0NL are zeroes, while their variances are related by 



d.S2(«(ri)«(r2)), (16) 



where both integrals arc along the same diffusion path. 
Because we have assumed nonlinearity to be weak, we can 
proceed by perturbation and replace the correlation func- 
tion of intensity fluctuations {SI{ri)SI{r2)) in Eq. lfTB|) 
by its value in the linear medium, i.e., by a sum of two 
contributions: the short-range one [Eq. ||2Jl] and the long- 
range one [Eq. (0J]. To perform integration in Eq. (|16|l . 
we replace Ar = |ri — in the expression jSJ for the 
short-range correlation function by As = |si — S2|j as- 
suming the wave path to be ballistic at distances shorter 
than £. In contrast, for Ar > £ the wave path is diffusive, 
and hence Ar in the expression for the long-range 
correlation function should be substituted by (As£)^/^. 
Performing then integrations in Eq. I|16|) yields 



(a</.?,l)^p(a<^^^. 



(17) 



with p the bifurcation parameter defined by Eq. H14|) with 
the first term, An'^{L/£)'^k£, resulting from the substitu- 
tion of the short-range correlation function ||2Jl in Eq. 
(|16|l and the second term, An'^{L/£)^, — from the sub- 
stitution of the long-range correlation function 10}. It 
follows from Eq. H17|) that at p > 1, the phase difference 
A(j) = A(/)L -f A0NL is dominated by A(/)nl and hence one 
may expect new phenomena (such as instability) to ap- 
pear. Because we have already shown that p = 1 is indeed 



the instability threshold, we adopt (A^^^) > (A^^) 



as 



the instability condition within our path-integral picture. 

The path-integral picture of wave propagation can be 
also applied to describe the dependence of instability 
threshold on the frequency of excitation. Indeed, assume 
that speckle pattern oscillates at some frequency that 
(we remind) is assumed to be much smaller than the in- 
verse mean free time c/£. At low frequencies Q ^ ^ 
we recover Eq. (|17|l because the oscillations are very slow 
and everything happens as if the speed of light diffusion 
through the random medium were infinite. If, in con- 
trast, fl > ^, situation changes. In Eq. I|16|) . the result 
of integration of the short-range correlation function Q 
remains the same as before since this correlation is only 
non-zero for Ar < £ and fi"^ is still much larger than the 
time needed for light to cover this distance. In contrast, 
integration of the long-range correlation function now 
yields a different result. Indeed, in this case correlation 
between intensities at ri and r2 can only be non-zero for 
Ar = |ri — r2| < {D^l~^y/^ < L due to the finite speed 
of light diffusion in random medium. Integration in Eq. 
(|16|l is therefore limited to pairs of points ri, r2 separated 
not more than by a distance (Z?r2~^)^/^. This gives 



L 



(A0!).(18) 



The condition of instability of speckle pattern with re- 
spect to excitations at frequency f2, (Ac/)^^) > (A'/'l)' 
can now be written as 



P 



> 



1 + k£/{L/£) 



[nTo)-^/^ + k£/{Ll£y 



(19) 



This coincides with Eq. IpUl for » > Tp ^ 

The path-integral picture of speckle pattern instabil- 
ity can help to understand how to apply the theoretical 
model developed in this paper to realistic experimental 
situations, when the disordered sample is characterized 
by three different size parameters L^, Ly, and (in the 
so-called 'slab' geometry, for example, = Ly ^ 
£). In this case, the only size parameter L of our theo- 
retical model should be understood as a size of the part 
of disordered sample explored by a typical diffuse path 
(and hence L « in the slab geometry). 



VII. DISCUSSION AND CONCLUSIONS 

Theoretical analysis performed in the present paper in- 
dicates that the multiple-scattering speckle pattern in a 
nonlinear random medium with intensity-dependent di- 
electric constant becomes unstable with respect to weak 
dynamic perturbations if the nonlinearity exceeds the 
threshold p ~ 1. This absolute instability threshold ap- 
pears to be independent of the nonlinearity response time 
TNL, whereas the frequency dependence of the threshold 
(and, in particular, the spectrum of frequencies subject 
to instability) is very sensitive to tnl- It is pertinent to 
note the extensive nature of the bifurcation parameter 
p. According to Eq. (|14|) . weak nonlinearity can be ef- 
ficiently compensated by a sufficiently large sample size 
L and p ^ 1 can be reached even at vanishing An, pro- 
vided that L is large enough. Another remarkable fea- 
ture of Eq. H14|l is the independence of p of the sign of 
An. The phenomenon of speckle pattern instability is 
therefore expected to develop in a similar way for both 
positive and negative nonlinear coefficients n2. This is 
not common for instabilities in nonlinear optical systems 
without disorder'^4,25,26,27,28,42,43,44,45,46,47 gj^^e the lat- 
ter instabilities are often related to the self-focusing ef- 
fect, arising at n2 > only. The instability discussed in 
the present paper is of different type and has nothing to 
do with self- focusing (which is negligible for Aii?k£ <Si 1 
assumed throughout the paper). This does not mean, 
however, that phenomena similar to those discussed in 
the present paper do not exist in homogeneous nonlin- 
ear media. In fact, it is easy to see that Eq. 1)15(1 de- 
scribes nothing but the self-phase modulations^ of the 
multiple-scattering speckle pattern. The development of 
self-phase modulation in a random medium, however, ap- 
pears to be rather different from that in the homogeneous 
case. Indeed, in a homogeneous medium of size L the 
nonlinear phase shift is deterministic and </)nl = kn2loL 
for a wave of intensity /q, whereas in the case of diffuse 
multiple scattering 0nl is random with the average value 
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((/■nl) ^ kn2lo{L'^/£) and variance (^nl) ~ (</'nl)^ ^ P- 
The threshold of the speckle pattern instability p ~ 1 can 
be therefore viewed as simply the point where the sample- 
to-samplc fluctuations of the nonlinear phase shift be- 
come of order unity. 

It is worth while to comment on the connection of the 
problem considered in the present paper to the recent re- 
sults on modulational instability (MI) of incoherent light 
beams in homogeneous nonlinear media— In this case, 
a 2D speckle pattern (speckled beam) with a controlled 
degree of spatial and temporal coherence is 'prepared' 
by sending a coherent laser beam on a rotating diffuser. 
The beam is then incident on a homogeneous nonlinear 
medium (inorganic photorefractive crystal). Pattern for- 
mation and 'optical turbulence' are observed if the inten- 
sity of the beam exceeds a specific threshold, determined 
by the degree of coherence of the beam. Although the 
study presented in this paper is also concerned with the 
instability of speckle patterns in nonlinear media, it dif- 
fers essentially from that of Refs.j^ and0 at least in the 
following important aspects: (a) light loses its spatial co- 
herence inside the medium, due to the multiple scattering 
on heterogeneities of refractive index, (b) the temporal 
coherence of incident light is perfect (in the case of in- 
coherent MI^S*^, coherence time of the incident beam is 
much shorter than the nonlinearity response time tnl), 
and (c) the incident light beam is completely destroyed 
by scattering after a distance ~ £ and light propagation 
is diffusive in the bulk of the sample. 

Let us now say a few words about practical conse- 
quences of our results and conditions of possible exper- 
imental observation of the instability phenomenon for 
light in random media. In a real experiment, infinitely 
weak dynamic perturbations of stationary speckle pat- 
tern are always present (due to, e.g., various types of 
noise). Our analysis predicts that the spectral compo- 
nents of these perturbations falling within the instability 
band (0, rimax) will be amplified due to the nonlinear- 
ity of the medium. Because our analysis is based on a 
linearized equation Q, we cannot predict further devel- 



opment of speckle dynamics with certainty, but our re- 
sults suggest the following scenario. In general, there are 
two possibilities for the spontaneous dynamics of speckle 
pattern beyond the threshold. First, it can be (quasi-) 
periodic, i.e. dominated by oscillations at some discrete 
frequency (or frequencies) Qi. Second, it can be chaotic 
with a continuous spectrum of intensity fluctuations be- 
tween some r^min and r^max- In nonlinear systems where 
dynamics is periodic, the frequency of oscillations fij is 
often favored already on the level of the linear stability 
analysis: it is often the frequency for which the maximum 
Lyapunov exponent A is reached. In our case, A exhibits 
monotonic decrease with (see Figs. 0] and 0), favoring 
no particular frequency. We have therefore no reason to 
expect the first scenario and the speckle pattern /(r, t) 
is likely to exhibit chaotic dynamics immediately beyond 
the instability threshold. The spectrum of chaotic fluctu- 
ations of /(r, t) will be continuous and comprised between 
^min — and f^max that we have discussed above. We 
note, however, that our analysis assumes diffusion regime 
of wave propagation {k£ :$> 1 and L/£^ 1) and is hkely 
to fail if the scattering is low-order {L < £) or if the 
regime of Anderson localization is approached {k£ ~ 1). 
In the latter cases, the first scenario of instability devel- 
opment can be a plausible option. 

Finally, experimental observation of instability of 
multiple-scattering speckle pattern will require, first of 
all, large enough sample size L and as low absorption as 
possible (if the macroscopic absorption length La < L, 
the effective size of the sample will be limited to La and 
one will not fully exploit the extensivity of instability 
threshold). In the absence of absorption (or for La ^ L), 
Att-nl 10^^ (realistic in nematic liquid crystals) and 
L/£ ^ 20 will suffice to get p ~ 1 and reach the insta- 
bility threshold. We note that in liquid crystals, optical 
nonlinearity is due to reorientation of molecules under 
the influence of electric field of electromagnetic wave and 
hence is essentially noninstantaneous. This emphasizes 
the importance of including the noninstantaneous nature 
of nonlinear response in our analysis. 
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